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Abstract 


The numerical simulation of steady planar two-dimensional, laminar flow of an 
incompressible fluid through an abruptly contracting channel using spectral domain 
decomposition methods is described. The key features of the method are the 
decomposition of the flow region into a number of rectangular subregions and 
spectral approximations which are pointwise C l wontinuous across subregion 
interfaces. Spectral approximations to the solution are obtained for Reynolds 
numbers in the range [0^ 500].” The size of the salient corner vortex decreases as the 
Reynolds number increases from 0 to around 45. As the Reynolds number is 
increased further the vortex grows slowly. A vortex is detected downstream of the 
contraction at a Reynolds number of around 175 that continues to grow as the 
Reynolds number is increased further. 

'Research was supported in part for the second author by the National 
Aeronautics and Space Administration under NASA Contract No. NAS1-1S605 while he 
was in residence at the Institute for Computer Applications in Science and 
Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 



1 . Introduction 


The numerical simulation of the flow in a constricted channel is considered 
using spectral domain decomposition techniques. The first numerical solutions to 
this problem were obtained by Dennis and Smith (1980) using a finite difference 
discretization of the stream function - vorticity formulation of the governing 
equations. They did not detect a downstream recirculation region caused by the 
flow separating at the corner despite using a very fine uniform grid. In a recent 
paper Hunt (1989) uses a non-uniform grid to ensure a dense distribution of grid 
points in the regions of the flow which need to be resolved. He places a locally 
fine mesh downstream of the constriction and around the re-entrant corner since 
this gives rise to a singularity in the vorticity. 

In this paper the governing equations are written in terms of the stream 
function. This means that mass is conserved identically. The flow region is 
divided into a number of rectangular subdomains. In each of these subdomains the 
stream function is approximated by a truncated double Chebyshev expansion. The 
expansion coefficients are determined by collocating the governing equation and 
appropriate interface continuity conditions between subdomains. The decomposition 
is such that the Chebyshev collocation points are distributed densely around the 
corner and near solid boundaries. This is important computationally in order to 
resolve the main features of the flow efficiently. 

In previous work, Karageorghis and Phillips (1989a) consider nonconforming 
subdomains because of their ease of implementation. Although this strategy works 
well for the Stokes problem, a lack of interface continuity appears for the Navier- 
Stokes problem (Karageorghis and Phillips (1989b)) for values of the Reynolds 
number around 200, eventually causing the method to break down completely. The 
spectral approximations obtained are not pointwise continuous across the subdomain 
interface. In this paper a collocation strategy is used which generates pointwise C 1 
continuous approximations across the interfaces. As a result no computational limit 
on the value of the Reynolds number has yet been encountered. Numerical 
solutions are obtained for values of the Reynolds number up to 500. 

The matrices arising from spectral discretizations are not sparse like their 
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finite difference and finite element counterparts. However, when used in 

conjunction with domain decomposition techniques these matrices possess a block 
tridiagonal structure which can be exploited when designing methods of solution. 
Here we use a subroutine from the NAG Library which solves almost block diagonal 
systems (Brankin and Gladwell (1990)). This method is as efficient as the 
capacitance matrix method in terms of cost and storage but is found to be more 
stable (Karageorghis and Phillips (1990)). 

Spectral approximations to the solution of this problem are obtained for 
Reynolds numbers in the range [0, 500]. A comparison with the work of Dennis and 
Smith (1980) shows good agreement between the sets of results in terms of stream 
function contours in the bulk of the flow and the description of the salient corner 
vortex. Although they did not find a downstream recirculation region there is a 
hint of its existence on their finest glides for Reynolds numbers above 1000. They 
probably failed to detect this behaviour in their numerical simulations because the 
grid employed was not fine enough in this region. In our calculations separation 
first appears at around Re - 175 which is slightly earlier than Hunt (1989) predicts. 
The separation length and the strength of the vortex increase as Re increases, the 
length being roughly proportional to Re. The spectral collocation method predicts 
longer separation lengths than Hunt (1989). On the finest grids used the spectral 
collocation method employs about a third of the number of degrees of freedom as 
Hunt (1989) and a fifteenth of the number of degrees of freedom as Dennis and 
Smith (1980). 

The Governing Equations 

We consider the steady laminar flow of an incompressible fluid through an 
abruptly contracting channel with walls at y - ±1 for x < o, y = ± i for x > 0 
and 1/2 ^ |y| <; 1 for x = 0. Upstream we impose parabolic Poiseuille flow and 
we suppose that the flow is parabolic again far enough downstream. Since the flow 
is symmetric about y = 0 it is only necessary to seek a solution for y 0. 

In terms of non-dimensional variables the incompressible Navier-Stokes 
equations are 


(v . 7)v Vp + (Re)" 1 7 2 v 


( 2 . 1 ) 
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where v = (u, v) is the velocity vector, p is the pressure and Re is the Reynolds 
number. The introduction of a stream function, #(x, y), defined by 


u 


dll) 

3y 



(2.3) 


means that the continuity equation (2.2) is satisfied identically and (2.1) becomes, 
after elimination of the pressure, 


— Re 


g £ - I 


3y 


(V 2 0) 




The boundary conditions are 


ip = l, |y=0 on y = 1, x ^ 0 and Y ~ \ , x ;> 0, 

ip = l, || = 0 on x = 0, i £ y £ 1 , 

ip = 0, =0 on y = 0. 

3y 


(2.4) 


(2.5) 


( 2 . 6 ) 


(2.7) 


The conditions on y = 0 preserve the symmetry of the flow while the other 
conditions represent the no-slip conditions. 

The governing equation (2.4) is a nonlinear partial differential equation for 
the stream function. It is solved in an iterative fashion using a Newton-type 
method to linearize it (Phillips (1984)). It is the linearized partial differential 
equation which is solved at each Newton step using a spectral collocation method. 
Let us write (2.4) in operator form, 

L«0 - 0 (2.8) 

where L is a nonlinear operator. Suppose that is some approximation to the 
solution of (2.8). To obtain an improved approximation we replace L by its 
linearization about and then solve the problem 

L '(#*)* = - L(^*) , 


(2.9) 
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where L'($) is the Frechet derivative of L at ip defined by 


L'W<(> = — Re 


jp (V-0) - I?? ^(V-0) 

3y 3x 3x 3y 


3x 




36 

3y 




( 2 . 10 ) 

The new approximation to the solution is thus ip* -r 4>* This completes a single 
Newton step and is repeated until convergence is reached. 


■L. Domain Decomposition and Spectral Approximation 

Since spectral methods are most easily applied to problems defined in 
rectangular regions, the natural domain decomposition of the flow region comprises 
three rectangular subdomains with common point (0, 1/2) as shown in Figure 1. Such 
a domain decomposition may be termed conforming since the sides of each rectangle 
are contiguous with either a boundary segment or a side of another rectangle. A 
further advantage of this decomposition is that it is possible to construct a 
collocation scheme which results in pointwise C 1 continuous approximations across 
the subdomain interfaces. 


The flow region is truncated upstream and downstream of the constriction at 
finite distances h 1 and h 2 from the origin, respectively. The distances h x and h 2 
need to be sufficiently large so that the flow is fully developed in the entry and 
exit sections. The domain truncation means that additional boundary conditions 
need to be imposed on entry and exit, namely that the normal derivative of ip 
vanishes on entry and exit. 


B C( 0,1) 


I 

D(0, 1/2) 

II 

III 


AC-hp 0) G(0, 0) F(h 2 , 0) 


Fig. 1 Three subdomain decomposition of flow region. 



- 6 - 


In each subdomain the stream function #(x, y) is approximated by V> k (x, y), 

where 

M k N k 

I& k (x, y) = G k (y) + a ™ n W " l(x) P " (y) ’ k = *’ n - IH - ai) 

m=M k n-N k 


and 


G ! (y) - G U (y) = \y(3-y-) , G HI (y) - G ! (2y); 

A 

,11 


Mo - N 0 


M 


0 


N 0 - M o" - No” - 4 • 


The polynomials (w£,(x)} and (P k (y)} are modified shifted Chebyshev polynomials 
which satisfy identically all the boundary conditions with the exception of the 
conditions along the vertical wall CD- For example, in region I we have 


wJn(x) = Tjn(x) + a l m t}(x) + AtJ(x) , 2 £ m £ M 1 , (3.2) 

where t{„(x) , 0 ^ m £ M 1 , are shifted Chebyshev polynomials on {— h^ 0] 
defined by 


tL(x) = 


t* r 2x + h n 
lm < iq J • 


and a{n and 0^ are given by 


oJ m = (-l) m m 2 , 0 { m = ( — l) m (m 2 - 1) , 2 £ m <; M 1 


Similarly, we can show that 

P* (y) = T [ n (y) + ajfj(y) + & n f l 0 iy) , 2 ^ n ^ N 1 , (3.3) 


where t\ (y) , 0 ^ n ^ , are the shifted polynomials on [1/2, 1] defined by 


fl(y) = 


T„ (4y-3) 
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and ctj , ]i\ are given by 

a 1 ,, = — n : — 1, 0l = n 2 — 1 . 

The basis functions in regions II and III are defined in a similar fashion. 


4. Collocation Strategy 

The coefficients {a^, n ) , k — « I, II, III must be determined for the spectral 

approximations (3.1) to be defined everywhere in the truncated flow region. 

Suppose that at the beginning of a Newton step we have approximations a** to the 

1c 

expansion coefficients of ip (x,y). The numerical solution of the linearized equation 
(2.9) determines the expansion coefficients (tfa)™ of <p y , the correction to lp v , in 
region k for k =1, II, III. We describe in some detail the process of calculating the 
coefficients (5a)Sm Once these are known, the updated approximation to the stream 
function is just ip* k + <p k , which can be found by simply adding together a** and 
($a)S, n for each value of m, n and k. 


The linearized partial differential equation (2.9) is collocated in each 
subdomain and the approximations in the three subregions are patched by imposing 
the correct order of weak continuity across the subregion interfaces. The 
collocation points in a subdomain are chosen to be those points at which the 
Chebyshev polynomial of highest degree used in the representation in that 
subdomain attains its extreme values. This choice gives rise to optimal 

approximation properties of smooth functions. It can also be shown that when 
Gaussian quadrature rules based on these points are used to evaluate the integrals 
appearing in Galerkin formulations of certain differential equations, then the 
resulting equations are equivalent to those determined by collocating the 
differential equation at the same set of points. Thus, in region I, for example, the 
collocation points are given by 


x 


I 

i 


h^Xj — 1) 
2 



y j + 3 


4 


(4.1) 


where 
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COS (~~j) , 0 <; i ^ M* : 

M 1 

cos (^) , 0 ^ j ^ N 1 . 

V 


Boundary conditions 

Due to the choice of modified Chebyshev polynomials as trial functions in 
the expansions (3.1) the boundary conditions are automatically satisfied except along 
x = 0 (1/2 <; y <£ 1) in region I. Along this portion of the boundary there are 
N^+l collocation points. We deduce from (3.1) that 0 and 0 X are polynomials of 
degree in y along x = 0 (1/2 <£ y 1) each depending on — 1 degrees of 
freedom. Therefore collocation of the boundary conditions at the — 1 distinct 
points (0, y*j), j = 0, 1, N 1 — 2, ensures that the no-slip boundary conditions 
are satisfied identically along this part of the boundary. 


Interf ace continuity conditions 

We impose C 3 continuity of the stream function across the interfaces 
between the subdomains. These conditions are collocated in such a way so as to 
yield approximations that are pointwise C 1 continuous across the interfaces, but 
whose second and third normal derivatives are continuous only at the interface 
collocation points. If the initial approximation to the solution of the problem 
satisfies these conditions then these continuity conditions are imposed on 0 at each 
Newton step. 

Let us examine the interface y — l/2( — h x x ^ 0) between subregions I 
and II, for example. Along y = 1/2, 0* and 0^ are polynomials of degree each 
possessing — 1 degrees of freedom. We collocate at a sufficient number of 
points on the interface to ensure that 0^ — 0^ is identically zero. This is 
achieved if we impose 

0 J (x | , i) — 0 II (x[ , i) = 0 , i - 2, 3, M 1 — 2, (4.2) 


and, in addition, 


0 n (O, i) = 


3 ^ 

3x 


(0, i) = 


(4.3) 


0 , 
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Similarly, continuity of <t>y is obtained by imposing 


II 


% <*! ■ b - • b - 


3y 


i = 2, 3, M 1 - 


(4.4) 


and 


3 ^ 

9y 


(0, i) = 


.LII 

O 0 

3x3y 


= 0 


(4.5) 


The conditions which represent the continuity of the second and third normal 
derivatives of 0 are collocated at the M* — 3 points (xj , 1/2) , i = 2, 3, 

M* — 2. Thus, these derivatives are not pointwise continuous across the interface. 
Moffatt (1964) shows that the leading term in the singular expansion of the Stokes 
solution about the re-entrant corner (0, 1/2) is O(r^) where X =* 1.5445 and so it 
would be inconsistent to impose pointwise continuity of these higher derivatives 
across the interface. 

The same collocation strategy is applied across the interface x =» 0 
(0 ^ y £ 1/2) between subregions II and III. As a result the correction 0 and its 
normal derivative are pointwise continuous across the interfaces. The updated 
approximations to the stream function also possess the same order of pointwise 
continuity. 


Linearized differential equation 

The linearized differential equation (2.9) is collocated in each subregion at all 
points on the collocation grid with the exception of those on or one in from each 
subregion boundary, i.e., 

(x k , y k ) , i=2, M k - 2, j = 2, .... N k - 2, 
for k - I, II or III. 


When the spectral collocation equations resulting from the boundary 
conditions, interface continuity conditions and differential equation are added 
together, they yield a total of KM^ — 1) (N* — 1) + (M^ — 1) (N^ — 1) + 
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(M^ — 1)(N^— 3)] linear equations which is equal to the number of unknown 
expansion coefficients (5a)^ n for k = I, II, III. Therefore, provided the coefficient 
matrix is non-singular, this system of equations possesses a unique solution. 


Direct Method of solution 

The global system for the expansion coefficients is 

Gz = r 

where 


A 

0 

0 


(5a) 1 

B 

0 



0 

c 

1 

[ 

0 

» z = 

(5a) 11 , 

0 

c 

f 

) 



0 

0 

E 


(5a) 111 

_ 


. 


. 


(5.1) 



(5.2) 


The vectors (5a)^, (5a)^ and (5a)^ contain the unknown coefficients in the 
expansions in subregions I, II and III, respectively. The first, third and fifth blocks 
of rows in (5.2) correspond to the collocation of the linearized differential equation 
and boundary conditions not satisfied by the approximations in subregions T, II, and 
III, respectively. The second and fourth blocks of rows correspond to the 


Hi II I 
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collocation of the interface continuity conditions between subregions I and II, and II 
and III, respectively. 

The system (5.1) is solved using an application of a standard production code 
designed for the solution of almost block diagonal systems (Brankin and Gladwell 
(1990)). In a recent survey of direct methods for solving systems of equations 
arising from spectral domain decomposition methods, Karageorghis and Phillips (1990) 
found this solver to be superior to other techniques with regard to cost, stability 
and storage. The code uses a modified column elimination procedure with alternate 
row and column pivoting based on an algorithm originally described in Varah (1976) 
and Diaz, Fairweather and Keast (1983) and is intended to solve systems of the form 
shown in Figure 2. These systems comprise rectangular blocks along the diagonal 
and are such that no three successive blocks have columns in common. 



Fig. 2 Almost block diagonal form. 

The global spectral collocation matrix G in (5.1) may be written in almost 
block diagonal form in the obvious way: 
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i 1 

i 

»— * 

o 

o 

1 

P 2 

1 q 2 

Q 1 

i R 2 



i 

0 

0 

R 2 



! r 3 

s 2 


T 

f -" 1 

1 

O 

j 

i 0 

T 2 


The matrix G has five non-zero blocks, namely, 



P 

Q 


1 


Q 2 

Ri 


[ R 2] 



and [T 2 ] . 


(5.3) 


Unfortunately the form (5.3) is not of the structure required by the almost block 
diagonal solver due to too much overlap between blocks 2 and 3, and 3 and 4. 
However, the transpose of G is of the required form, with three blocks, namely, 


[p T I Qj , Qj I R T l S J ] and 



One may, therefore, decompose the transpose of the global matrix using the 
existing NAG subroutine F01LHF and subsequently solve using the transpose of the 
decomposed form of G , say G 1 , the system 



(5.4) 


with the NAG subroutine F04LHF. Further details of the implementation may be 
found in Karageorghis and Phillips (1990). 
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6^ Numerical Results 

The results of numerical calculations of laminar flow through a 2:1 

constricted channel are presented and comparisons made with the work of Dennis and 
Smith (1980) and Hunt (1989). Dennis and Smith (1980) use a finite difference 
discretization of the stream function-vorticity formulation of the Navier-Stokes 
equations on a uniform grid. An upwind differencing scheme of Dennis and Hudson 
(1978) is used to approximate the vorticity transport equation. This is essential for 
reasonably high values of Re in order to maintain the diagonal dominance in the 
approximating sets of difference equations. The loss of diagonal dominance presents 
difficulties in the convergence of iterative solution techniques. If upwinding is not 
used then diagonal dominance can only be established either by sufficiently 

decreasing the mesh size, which may be prohibitively expensive, or by using a non- 
uniform grid. Hunt (1989) also uses a finite difference discretization of the stream 
function-vorticity formulation but on a non-uniform grid. The vorticity unknowns 
are eliminated to give a system solely in terms of unknown values of the stream 
function. Hunt (1989) investigates the application of artificial viscosity or 
upwinding but contrary to expectations finds that the scheme with upwinding fails 
to converge for Re > 500 even though one would expect the addition of a Viscous’ 
like term to have a stabilizing effect. Further for a given value of Re convergence 
becomes increasingly more difficult on successively finer grids. The opposing 
conclusions on the application of upwind differencing must be due to the different 
nature of the grids used in these studies. The ability of these two methods and 
the spectral algorithm to describe the main features of the flow is examined. 

The solution of the Stokes problem (Re = 0) is used as the initial 

approximation to the solution of the Navier-Stokes problem for 0 < Re ^ 150. 
However, for Re <; 150 the Newton process is robust enough to converge from an 
initial approximation that is zero everywhere. For Re > 150 continuation in Re is 
required for convergence in increments of Re of 50. This means, for example, that 
the converged numerical solution for Re = 150 is used as the initial approximation 
when calculating the numerical solution for Re = 200. 

The Newton process is terminated when the maximum residual is less than 
10- 7 in magnitude. This invariably occurs after six Newton iterations. The 
expansion coefficients are also checked for convergence. The last two iterations 
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are found not to affect the first eight significant digits of all the expansion 
coefficients. 

Numerical experiments are performed on a number of collocation grids. This 
is necessary in order to ensure that the approximations obtained are not strongly 
grid dependent. The following discretization parameters define the collocation 
grids used in this work and the total number of degrees of freedom: 

(1) M 1 = M 11 = 24, N 1 = 20, N 11 = N 111 = 16, M 111 = 16 (977 degrees of 
f reedom); 

(2) M 1 = M 11 = 24, N 1 = N 11 = N 111 = 20, M 111 = 24 (1265 degrees of 
freedom); 

(3) M 1 = M 11 - 24, N 1 - N 11 = N I]I = 20, M 111 = 32 (1401 degrees of 
freedom); 

(4) M 1 = M 11 = 24, N 1 = N 11 = N 111 = 20, M 111 = 40 (1537 degrees of 
freedom). 


For Re <; 100 good agreement is obtained on all four grids. For 100<;Re^300 
good agreement is reached on grids (2) — (4). Thereafter, for Re ^ 300, grids (3) 
and (4) produced almost identical results. 

The upstream and downstream truncation lengths are chosen to be 1.0 and 3.5, 
respectively. It is necessary to experiment with the number of degrees of freedom 
in subregion III since for an insufficient number the flow is not adequately resolved 
particularly near the downstream recirculation region. Dennis and Smith (1980) 
admit that they do not resolve this feature of the flow and that the use of very 
fine grids is necessary to recover the true behaviour. However, from their 
numerical calculations they are able to imply that the fluid separates downstream of 
the constriction. 

The behaviour of the vorticity along the downstream channel wall (y = j, x 

2> 0) indicates whether the domain has been truncated far enough downstream. As 

2 1 

x - oo the vorticity f = -V # -> 12. If f(x, ^) settles down to this value then the 
truncation length is adequate. The value of f(x, ±) is plotted in Figs. 3(a) — (c) 
for Re = 0, 100, 500, respectively. These figures suggest that the exit length is 
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suitable for the range of Reynolds numbers considered here. Dennis and Smith 
(1980) use downstream truncation lengths of 2 and 2.5 and obtain excellent agreement 
in their results which demonstrates that the numerical solutions obtained with a 
truncation length of 2 are satisfactory. Hunt (1989) uses a transformation of the 
independent variables to set the downstream boundary at x — 1000. Instead of 
experimenting with truncation lengths he needs to determine appropriate mapping 
parameters which appear in the transformation. 

Contours of the stream function are presented in Fig. 4(a) — (i) for Re = 0, 
10, 50, 100, 150, 200, 300, 400 and 500, respectively, in the region — 1 <; x <; 1, 
0 ^ y 1* The vortex in the salient corner diminishes in size from Re-0 until 
around Re = 45 and then increases slowly as Re is increased further. Let L v 
denote the distance between the point where the separation line meets the top of 
the channel and the salient corner. The value of L v is recorded in Table I for 
different methods and for a range of values of Re. The results of Dennis and 
Smith (1980) are obtained by means of two successive h 2 — extrapolation operations 
on information calculated on grids of mesh lengths 1/10, 1/20 and 1/40. Hunt (1989) 
uses a transformed grid with 48 x 128 points. The values of L v in columns (b) and 
(c) agree to within 5%. For the range of values of Re considered here the scheme 
of Hunt (1989) that uses artificial viscosity gives results that are closer to those in 
column (b) than the scheme without artificial viscosity. 

Close-up views, — - <; x <; 0, / <; y <; 1, of the salient corner are given in 
Fig. 5(a) — (e) for Re = 10, 50, 100, 150 and 200, respectively. The first closed 
streamline corresponds to a contour height of 1 + 10“ 5 and the remaining ones differ 
by multiples of 10' 5 . The interesting feature in these plots is the second vortex 
appearing close to the corner. In the Stokes case Moffatt (1964) predicts an infinite 
sequence of eddies running into the corner. The size of this second vortex grows 
moderately as Re is increased. Dennis and Smith (1980) only just detect the second 
vortex at values of Re of at least a thousand and then only on extremely fine 
meshes of size 1/80. 

A small recirculation region downstream of the constriction is observed in 
Fig. 4 (d) for Re = 100. This region suddenly grows when Re is increased to 200. 
In fact when the stream function is calculated at intermediate values of Re this 
downstream recirculation region grows suddenly at a value of Re around 175. A 
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magnification of this region is shown in Fig. 6 (a) - (d) for Re - 200, 300, 400 and 
500, respectively. In these figures the stream function is contoured in the 
rectangle 0 <. x <, 1, 1 <[ y <. v These figures give a good description of the way 
in which the vortex develops as Re is increased. In Fig. 6 the fluid is seen to 
separate at the corner and not slightly to the right as Hunt (1989) observes. 
Further the length of the downstream recirculation regions are significantly longer 
than those calculated by Hunt (1989) as can be seen in Table II. Dennis and Smith 
(1980) do not detect this region although they say that its existence is implied as 
the grid is refined. For Re — 500 they predict that separation occurs at a point in 
the interval 0 < x < 0.3 with reattachment at a point beyond x - 1.2. For this 
value of Re their prediction of the reattachment point is closer to our calculated 
value than the value obtained by Hunt (1989). 

7. Concluding Remarks 

The steady planar two-dimensional laminar flow of an incompressible fluid 
through an abruptly contracting channel is considered for moderate values of the 
Reynolds number. The governing equation for the stream function is linearized 
using Newton’s method and solved using spectral domain decomposition techniques. 
A conforming domain decomposition is used which in conjunction with a carefully 
constructed collocation strategy ensures that the resulting spectral approximations 
are globally - continuous. In addition, the spectral approximations are C°° 
except along subdomain interfaces. An efficient direct method for solving the 
spectral collocation equations at each stage of the Newton process is described. 

At most six Newton iterations are required for convergence. For 0 < Re <[ 
150 a converged numerical solution is obtained from an initial approximation that is 
zero everywhere. Continuation in Re in increments of 50 is used for 
150 < Re <[500. The vortex in the salient corner decreases in size from Re - 0 to 
around Re - 45 and then grows slowly as Re is increased further. A small 
recirculation region just downstream of the constriction appears at Re - 100. This 
region suddenly grows as Re is increased to a value around 175. This recirculation 
region extends further downstream as Re is increased. 

There is qualitative and quantitative agreement with Dennis and Smith (1980) 
in the bulk of the flow and the description of the salient corner vortex. Their 
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numerical solutions do not predict a downstream recirculation region although they 
say there is a hint of its existence as the grid is refined. Our calculations predict 
that the fluid separates at the reentrant corner not just to the right of it as does 
Hunt (1989) and that reattachment occurs further downstream than in his 
simulations. 

This paper demonstrates that spectral domain decomposition techniques are 
capable of solving complex flow situations and resolving the main features of the 
flow. This is accomplished in an efficient manner using far fewer degrees of 
freedom than other methods of discretization. 


ACKNOWLEDGEMENT 

The authors wish to thank Dr A. A. Minzoni for many enlightening 
discussions. 



- 18 - 


REFERENCES 


BRANKIN, R.W. & GLADWELL, I. 1990 Codes for almost block diagonal systems. 

To appear in Comp. Atath. Applies. 

DENNIS, S.C.R. & HUDSON, J.D. 1978 A difference method for solving the 
Navier-Stokes equations. Proceedings of the First Conference on 
Numerical Methods in Laminar and Turbulent Flow, Pentech Press, London. 

DENNIS, S.C.R. & SMITH, F.T. 1980 Steady flow through a channel contraction 
with a symmetrical constriction in the form of a step. Proc. Roy. Soc. 

London A 372, 393-414. 

DIAZ, J.C., FAIRWEATHER, G. & KEAST, P. 1983 FORTRAN packages for solving 
certain almost block diagonal linear systems by modified alternate row 
and column eliminations. ACM Trans. Math. Softw. 8, 358-375. 

HUNT, R. 1989 The numerical solution of the laminar flow in a constricted 
channel at moderately high Reynolds number using Newton iteration. 

Submitted to Int. J. Numer. Aleth Fluids. 

KARAGEORGHIS, A. & PHILLIPS, T.N. 1989a Spectral collocation methods for 
Stokes flow in contraction geometries and unbounded domains. J. Comput. 
Phys. 80, 314-330. 

KARAGEORGHIS, A. & PHILLIPS, T.N. 1989b Chebyshev spectral collocation 

methods for laminar flow through a channel contraction. J. Comput. Phys. 84, 
114-133. 

KARAGEORGHIS, A. & PHILLIPS, T.N. 1990 On efficient direct methods for 
conforming spectral domain decomposition techniques. Submitted to J. 

Comput. Appl. Math. 

MOFFATT, H.K. 1964 Viscous and resistive eddies near a sharp corner. J. Fluid 
Mech. 18, 1-18. 

NUMERICAL ALGORITHMS GROUP LIBRARY Mark 13, NAG (UK) Ltd., Banbury 
Road, Oxford, UK. 

PHILLIPS, T.N. 1984 Natural convection in an enclosed cavity. J. Comput. 

Phys. 54, 365-381. 

PHILLIPS, T.N. & KARAGEORGHIS, A. 1989 A conforming spectral collocation 
strategy for Stokes flow through a channel contraction. ICASE Report 
No. 89-60. Submitted to Appl. Numer. Aiath. 



- 19 - 


VARAH, J.M. 1976 Alternate row and column elimination for solving certain 
linear systems. 57. -1.1/ J. Numer. Anal. 13, 71-75. 



- 20 - 


TABLE I 


Re 

(a) 

(b) 

(c) 

(d) 

(e) 

0 

0.285 

0.290 

0.284 

- 

- 

10 

0.150 

0.148 

0.155 

- 

- 

25 

0.129 

0.125 

- 

- 

- 

50 

0.129 

0.123 

0.129 

- 

- 

75 

0.135 

0.135 

- 

• 

- 

100 

0.143 

0.140 

0.144 

- 

- 

125 

0.154 

0.150 

- 

0.164 

0.168 

150 

0.160 

0.155 

- 

- 

- 

200 

0.183 

0.183 

- 

- 

- 

250 

0.210 

0.205 

- 

0.209 

0.227 

300 

0.223 

0.223 

- 

- 

- 

350 

0.235 

0.233 

- 

- 

- 

400 

0.244 

0.244 

- 

- 

- 

450 

0.251 

0.255 

- 

- 

- 

500 

0.260 

0.265 

0.266 

0.260 

0.308 


The length of the upstream vortex as a function of Re for (a) spectral collocation 
method on grid (3), (b) spectral collocation method on grid (4), (c) extrapolated finite 
difference scheme of Dennis and Smith (1980), (d) finite difference scheme of Hunt 
(1989) with artificial viscosity, (e) finite difference scheme of Hunt (1989) without 
artificial viscosity. 
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TABLE II 


Re 

(a) 

(b) 

(c) 

150 

0.088 

- 

- 

200 

0.363 

- 

- 

250 

0.500 

0.106 

0.096 

300 

0.538 

- 

- 

350 

0.700 

- 

- 

400 

0.738 

- 

- 

450 

0.924 

- 

- 

500 

0.995 

0.406 

0.406 


The length of the downstream vortex as a function of Re for (a) spectral 
collocation method on grid (4), (b) finite difference scheme of Hunt (1989) with 
artificial viscosity, (c) finite difference scheme of Hunt (1989) without artificial 
viscosity. 



CAPTIONS FOR FIGURES 


Fig. 3. Wall vortcity f(x,l/2) for x > 0, for various Reynolds numbers Re. 

Fig. 4. Streamfunction contours in -1 <> x ^ 1, 0 ^ y ^ 1, for various Reynolds 
numbers Re. 

Fig. 5. Close-up views of the salient corner streamlines on -1/8 <; x <; 0, 

7/8 y <; 1 for various Reynolds numbers Re. 

2 1 

Fig. 6. Close-up views of the streamlines in the region 0^x<;l,^<;y<;^ 
for various Reynolds numbers Re. 
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Fig. 3(c) 
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Fig. 5 (b) 
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Fig. 5(e) 
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